{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Kennard-Stone 采样方法在 Python/C 下的高效实现 $O(N^2)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> 创建时间：2021-01-29"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在这篇文档中，我们会非常简单地回顾 Kennard-Stone 采样方法的原理，并集中精力讨论其算法实现方式。\n",
    "\n",
    "标题的 $O(N^2)$ 的 $N$ 表示的是样本数量 $n_\\mathrm{sample}$。严格的计算复杂度是 $O(n_\\mathrm{sample}^2 n_\\mathrm{feature} + n_\\mathrm{sample} n_\\mathrm{result})$，其中 $n_\\mathrm{feature}$ 为特征向量长度且认为是常数值，$n_\\mathrm{result}$ 为采样数。\n",
    "\n",
    "Kennard-Stone Algorithm 的原始文献是 Kennard, Stone [^Kennard-Stone.T.1969]。\n",
    "\n",
    "Kennard-Stone 采样的重要性在于，对于绝大多数机器学习方法而言，选择一个较好的采样方式构建训练集，会一定程度地增强泛化能力。该采样方法在小样本数据集中较常使用；大样本数据集大多采用随机分割方法。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这篇文档的特色有四处：\n",
    "\n",
    "- 通过 Python 中的 numpy 函数，可以并行、低额外内存消耗、且快速地构建欧式距离 (Euclidean Distance) 矩阵；\n",
    "\n",
    "- 通过活用 numpy 的函数，可以高效串行地实现 Kennard-Stone 采样方法，且相信比目前大多数开源方案快很多；\n",
    "\n",
    "- 通过 C/OpenMP 与 ctype binding，可以并行且进一步压榨采样效率。\n",
    "\n",
    "- 通过高效的程序实现，对于 Intel-i5 7300HQ 的 4 核个人电脑，30000 个样本 (特征向量长度为 100) 的采样时间最快仅仅需要不到 10 秒。如果有更多的内存，这几乎意味着十万 (100k) 级别的数据可以迅速以 Kennard-Stone 方法采样。\n",
    "\n",
    "需要注意，C 语言的程序需要在 gcc 的编译器上实现。同时，计算过程采用了单浮点，因此可能会在超级大样本中产生结果的差异。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:05.850456Z",
     "iopub.status.busy": "2021-01-31T13:14:05.850050Z",
     "iopub.status.idle": "2021-01-31T13:14:05.913935Z",
     "shell.execute_reply": "2021-01-31T13:14:05.913436Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.ma as ma\n",
    "\n",
    "np.set_printoptions(precision=6, linewidth=120, suppress=True)\n",
    "np.random.seed(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Python 实现"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 随机数据集"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "首先我们创造一个随机数据集，该数据集包含 5000 个样本。\n",
    "\n",
    "- `n_sample` $n_\\mathrm{sample}$ 为样本数。现在定为 5000。\n",
    "\n",
    "- `n_feature` $n_\\mathrm{feature}$ 为特征数 (描述样本性质的数据)。现在定为 100。\n",
    "\n",
    "- `X` $\\mathbf{X}$ 为完整的数据集。它是一个 $(n_\\mathrm{sample}, n_\\mathrm{feature})$ 大小的矩阵。\n",
    "\n",
    "    - 单个数据写为向量 $\\boldsymbol{x}_{i}$，矩阵元写为 $x_{ia}$。其中，$i, j, k, \\cdots$ 为样本角标，$a$ 为特征角标。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:05.916823Z",
     "iopub.status.busy": "2021-01-31T13:14:05.916432Z",
     "iopub.status.idle": "2021-01-31T13:14:05.931241Z",
     "shell.execute_reply": "2021-01-31T13:14:05.930846Z"
    }
   },
   "outputs": [],
   "source": [
    "n_sample  = 5000\n",
    "n_feature = 100\n",
    "X = np.random.randn(n_sample, n_feature)\n",
    "X *= 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由于这份文档的目标是说明代码的实现方式与效率，因此不打算具体地用小数据作结果演示，而是用使用中等的数据集。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 距离矩阵的生成 (欧氏距离)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "两个样本之间的欧氏距离 (Euclidean Distance) 定义为\n",
    "\n",
    "$$\n",
    "d_{ij} = \\Vert \\boldsymbol{x}_i - \\boldsymbol{x}_j \\Vert_2 = \\sqrt{\\sum_a (x_{ia} - x_{ja})^2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "因此，距离矩阵 $\\mathbf{D}$ 会是一个维度为 $(n_\\mathrm{sample}, n_\\mathrm{sample})$。\n",
    "\n",
    "但如果直接用上式计算欧式距离，计算代价会稍大。因此，一个取巧的方法是利用矩阵相乘的结果。\n",
    "\n",
    "定义中间变量的矩阵 $\\mathbf{Y} = \\mathbf{X} \\mathbf{X}^\\dagger$，并定义向量 $\\boldsymbol{t} = \\mathrm{diag} (\\mathbf{Y})$ 即 $\\mathbf{Y}$ 的对角线，或写成\n",
    "\n",
    "$$\n",
    "y_{ij} = \\sum_{a} x_{ia} x_{ja}, \\quad t_i = y_{ii}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "那么距离矩阵的平方就可以写为\n",
    "\n",
    "$$\n",
    "d_{ij}^2 = \\sum_a (x_{ia} - x_{ja})^2 = \\sum_a x_{ia}^2 - \\sum_a x_{ia} x_{ja} + \\sum_a x_{ia}^2 = t_i - 2 y_{ij} + t_j\n",
    "$$\n",
    "\n",
    "这就可以很大程度上简化计算量，并且不会引入很大的内存损耗。对于 numpy，由于 Boardcasting 机制，内存消耗可以是严格的 $n_\\mathrm{sample}^2 + n_\\mathrm{sample}$。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "计算复杂度是 $O(n_\\mathrm{sample}^2 n_\\mathrm{feature})$，这是因为两 $(n_\\mathrm{sample}, n_\\mathrm{sample})$ 维度的矩阵相乘复杂度便是如此。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在我们考察的实现过程：\n",
    "\n",
    "- `dist` 距离矩阵 $\\mathbf{D}$ 或写成 $d_{ij}$，维度 $(n_\\mathrm{sample}, n_\\mathrm{sample})$；它是这一步的输出。\n",
    "\n",
    "- `t` 中间变量矩阵 $\\boldsymbol{t} = \\mathrm{diag} (\\mathbf{Y})$ 或写成 $t_i$，维度 $(n_\\mathrm{sample}, )$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:05.934988Z",
     "iopub.status.busy": "2021-01-31T13:14:05.934532Z",
     "iopub.status.idle": "2021-01-31T13:14:05.936949Z",
     "shell.execute_reply": "2021-01-31T13:14:05.936537Z"
    }
   },
   "outputs": [],
   "source": [
    "def get_dist(X):\n",
    "    dist = X @ X.T              # - 1\n",
    "    t = dist.diagonal().copy()  # - 2\n",
    "    dist *= -2                  # - 3\n",
    "    dist += t[:, None]          # - 4\n",
    "    dist += t[None, :]          # - 5\n",
    "    return np.sqrt(dist)        # - 6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "其具体算法过程是\n",
    "\n",
    "1. 生成 $y_{ij} = \\sum_a x_{ia} x_{ja}$，并将 $y_ij$ 储存在 `dist` 中；\n",
    "\n",
    "2. 生成 $t_i = y_{ii}$，并将 $t_i$ 储存在 `t` 中；\n",
    "\n",
    "3. 生成 $- 2 y_{ij}$，覆盖在 `dist` 中；\n",
    "\n",
    "4. 生成 $- 2 y_{ij} + t_i$，覆盖在 `dist` 中；\n",
    "\n",
    "5. 生成 $- 2 y_{ij} + t_i + t_j$，覆盖在 `dist` 中；\n",
    "\n",
    "6. 生成 $\\sqrt{- 2 y_{ij} + t_i + t_j}$，覆盖在 `dist` 中。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "那么当前样本所对应的距离矩阵就通过下述代码生成："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:05.940001Z",
     "iopub.status.busy": "2021-01-31T13:14:05.939605Z",
     "iopub.status.idle": "2021-01-31T13:14:06.165926Z",
     "shell.execute_reply": "2021-01-31T13:14:06.166445Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 592 ms, sys: 61.2 ms, total: 654 ms\n",
      "Wall time: 223 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "dist = get_dist(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最后指出，这部分代码显然应当可通过 C 语言加速，但由于 numpy 的底层计算实现已经利用了 C 语言与经过优化的 Blas，因此我认为没有必要对距离矩阵的生成再使用 C 语言优化。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 默认的种子样本：距离最远的两个样本"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Kennard-Stone 采样需要基于一部分“种子样本”(Seed) 进行。令记号 $s_i$ 代表种子样本，$n_\\mathrm{seed}$ 代表种子数量。\n",
    "\n",
    "默认的种子是距离最远的两个样本，即 $n_\\mathrm{seed} = 2$，\n",
    "\n",
    "$$\n",
    "(s_0, s_1) = \\arg \\max_{(i, j)} d_{ij} \\text{, or equilvently, } d_{s_0 s_1} = d_{s_1 s_0} = \\max_{(i, j)} d_{ij}\n",
    "$$\n",
    "\n",
    "寻找最大的元素位置可以直接地通过 `numpy.argmax` 函数得到。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下面定义一些程序变量：\n",
    "\n",
    "- `max_indexes`：距离最远的两个样本指标 $\\arg \\max_{(i, j)} d_{ij}$，为 tuple 型变量；\n",
    "\n",
    "- `max_dist`：最远的距离 $d_\\mathrm{max} = \\max_{(i, j)} d_{ij}$；\n",
    "\n",
    "- `seed`：种子样本 $\\boldsymbol{s}$ 或写为 $s_i$，维度为 $(n_\\mathrm{seed}, )$；它原则上可以是任何非空样本集，但默认是 $(s_0, s_1) = \\arg \\max_{(i, j)} d_{ij}$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:06.169033Z",
     "iopub.status.busy": "2021-01-31T13:14:06.168261Z",
     "iopub.status.idle": "2021-01-31T13:14:06.199013Z",
     "shell.execute_reply": "2021-01-31T13:14:06.199545Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 78.7 ms, sys: 0 ns, total: 78.7 ms\n",
      "Wall time: 19.7 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(450, 4092)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "max_indexes = np.unravel_index(np.argmax(dist), dist.shape)\n",
    "max_dist = dist[max_indexes]\n",
    "seed = max_indexes\n",
    "seed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Kennard-Stone 采样"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "令 $r_0, r_1, \\cdots, r_{n - 1}$ 代表的是第 $0, 1, \\cdots, n-1$ 次采样的样本序号。现在需要采样第 $n$ 个样本。\n",
    "\n",
    "如果我们现在知道，对于任意**未采样**的第 $i$ 个样本，在第 $n$ 步时有最小距离向量\n",
    "\n",
    "$$\n",
    "{}^n m_i = \\min \\{ d_{r_0 i}, d_{r_1 i}, \\cdots, d_{r_{n-1} i} \\}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "那么第 $n$ 步被采样的样本应当是 (到已采样点最小距离能取到最大值的样本)\n",
    "\n",
    "$$\n",
    "r_n = \\underset{i \\not \\in \\{ r_0, r_1, \\cdots, r_{n-1} \\}}{\\text{arg max}} \\big\\{ \\min \\{ d_{r_0 i}, d_{r_1 i}, \\cdots, d_{r_{n-1} i} \\} \\big\\} = \\underset{i \\not \\in \\{ r_0, r_1, \\cdots, r_{n-1} \\}}{\\text{arg max}} {}^n m_i\n",
    "$$\n",
    "\n",
    "需要注意，上式的 $i$ 必须是在**未采样的样本**中得到的；但这个记号较为冗余，因此后文会省略。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "随后需要为下一步采样过程作准备。下一步需要使用到更新后的最小距离向量\n",
    "\n",
    "$$\n",
    "{}^{n+1} m_i = \\min \\{ d_{r_0 i}, d_{r_1 i}, \\cdots, d_{r_{n-1} i}, d_{r_n i} \\}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这个向量的构建对于单个数据 $i$ 来说，不需要 $O(n)$ 的计算复杂度，因为上式等价于\n",
    "\n",
    "$$\n",
    "{}^{n+1} m_i = \\min \\{ {}^{n} m_i, d_{r_n i} \\}\n",
    "$$\n",
    "\n",
    "随后就可以得到 $r_{n+1}$，以此类推。这是算法可以变得较快的实际原因。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Kennard-Stone 算法的简述就是这些。对于程序作下述说明。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**函数输入与输出**\n",
    "\n",
    "已经定义的量是\n",
    "\n",
    "- 输入量 `dist` 距离矩阵 $d_{ij}$，维度 $(n_\\mathrm{sample}, n_\\mathrm{sample})$\n",
    "\n",
    "- 输入量 `seed` 种子样本 $s_i$，维度 $(n_\\mathrm{seed}, )$。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "还需要定义的\n",
    "\n",
    "- 输入量 `n_result` $n_\\mathrm{result}$ 待采样数据数量；这里我们进行 4990 个样本采样；\n",
    "\n",
    "- 输出量 `result` $\\boldsymbol{r}$ 或 $r_i$ 为采样结果序号，维度 $(n_\\mathrm{result}, )$，并依照采样顺序记录结果；\n",
    "\n",
    "- 输出量 `v_dist` $\\boldsymbol{v}$ 或 $v_i$ 为采样距离，维度 $(n_\\mathrm{result}, )$：\n",
    "\n",
    "    $$\n",
    "    v_{i} = \\max_j {}^{i + 1} m_j = \\max_j \\big\\{ \\min \\{ d_{r_0 j}, d_{r_1 j}, \\cdots, d_{r_i j} \\} \\big\\}\n",
    "    $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**中间变量定义** (`Definition: Intermediate Variables`)\n",
    "\n",
    "- `n_seed` 种子样本的维度 $n_\\mathrm{seed}$；\n",
    "\n",
    "- `selected` 储存是否已被采样的信息向量，维度 $(n_\\mathrm{sample}, )$；\n",
    "\n",
    "- `min_vals` 最小值向量 ${}^n \\boldsymbol{m}$ 或 ${}^n m_i$，维度 $(n_\\mathrm{sample}, )$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**初始化过程** (`--- Initialization ---`)\n",
    "\n",
    "1. 首先将种子 $(s_0, s_1, \\cdots, s_{n_\\mathrm{seed}})$ 放到结果中，即令 $r_0 = s_0$, $r_1 = s_1$, $\\cdots$, $r_{n_\\mathrm{seed}} = s_{n_\\mathrm{seed}}$；\n",
    "\n",
    "2. 如果种子样本数量恰为 2，那么就认为选取了欧氏距离最远的两个点，并认为 $v_0 = d_{s_0 s_1} = \\max_{(i, j)} d_{ij}$；\n",
    "\n",
    "3. 初始化种子样本下的最小距离向量。注意到 $r_0 = s_0$，因此\n",
    "\n",
    "    $$\n",
    "    {}^1 m_i = \\min \\{ d_{r_0 i} \\} = d_{s_0 i}\n",
    "    $$\n",
    "\n",
    "4. 确定一个值 `upper_bound`，使得对于任意的 $n, i$，都有 ${}^n m_i \\leqslant \\mathtt{upper\\_bound}$。这是利用了 ${}^n m_i \\leqslant {}^n m_i \\leqslant \\max_i {}^n m_i = \\mathtt{upper\\_bound}$ 的性质而写的\n",
    "\n",
    "    - 这个变量存在与否不影响算法本身，它是出于 numpy.min 函数在含有 `where` 传参时，必须要指定一个 `initial` 参数而定；\n",
    "\n",
    "5. 对于所有未选中的样本，都使用下述方式迭代更新最小距离向量：\n",
    "\n",
    "    $$\n",
    "    {}^{n+1} m_i = \\min \\{ {}^{n} m_i, d_{r_n i} \\} = \\min \\{ {}^{n} m_i, d_{s_n i} \\}\n",
    "    $$\n",
    "    \n",
    "    注意如果 $i$ 已经在被采样样本中，就无需作更新了；通过 `selected` 确定是否要更新样本 $i$。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**循环过程** (`--- Loop argmax minimum ---`)\n",
    "\n",
    "循环的目的是确定第 $n$ `n` 个样本所在位置。$n$ 从种子数 $n_\\mathrm{seed}$ 开始，到采样数 $n_\\mathrm{result}$ 结束。\n",
    "\n",
    "1. 获得采样的样本位置 $r_n = \\arg \\max_i {}^n m_i$，并储存在临时变量 `sup_index` 中；\n",
    "\n",
    "2. 将样本位置放到向量 $\\boldsymbol{r}$ 中；\n",
    "\n",
    "3. 具体的数值 $\\max_i {}^n m_i$ 放到 $\\boldsymbol{v}$ 中；\n",
    "\n",
    "4. 告知 `selected` 样本 $i$ 已经被选中了；\n",
    "\n",
    "5. 更新最小距离向量：${}^{n+1} m_i = \\min \\{ {}^{n} m_i, d_{r_n i} \\}$。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "程序所需要的额外内存是 $2 n_\\mathrm{result} + 2 n_\\mathrm{sample}$；这相对于 $n_\\mathrm{sample}^2$ 的距离矩阵来说已经不是很重要了。\n",
    "\n",
    "程序的计算复杂度是 $O(n_\\mathrm{sample} n_\\mathrm{result})$。这是因为需要采样 $n_\\mathrm{result}$ 个样本；每次需要在 $n_\\mathrm{sample}$ 数量的最小距离向量中寻找最大值，以及更新最小距离向量，因此复杂度就是采样数与样本数的简单乘积。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:06.202571Z",
     "iopub.status.busy": "2021-01-31T13:14:06.201785Z",
     "iopub.status.idle": "2021-01-31T13:14:06.205408Z",
     "shell.execute_reply": "2021-01-31T13:14:06.205938Z"
    }
   },
   "outputs": [],
   "source": [
    "n_result = 4990\n",
    "result = np.zeros(n_result, dtype=int)\n",
    "v_dist = np.zeros(n_result, dtype=float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:06.208495Z",
     "iopub.status.busy": "2021-01-31T13:14:06.207697Z",
     "iopub.status.idle": "2021-01-31T13:14:06.214339Z",
     "shell.execute_reply": "2021-01-31T13:14:06.214886Z"
    }
   },
   "outputs": [],
   "source": [
    "def ks_sampling_core(dist, seed, n_result):\n",
    "    # Definition: Output Variables\n",
    "    result = np.zeros(n_result, dtype=int)\n",
    "    v_dist = np.zeros(n_result, dtype=float)\n",
    "    # Definition: Intermediate Variables\n",
    "    n_seed = len(seed)\n",
    "    selected = np.zeros(n_sample, dtype=bool)\n",
    "    min_vals = np.zeros(n_sample, dtype=float)\n",
    "    # --- Initialization ---\n",
    "    result[:n_seed] = seed                   # - 1\n",
    "    if n_seed == 2:\n",
    "        v_dist[0] = dist[seed[0], seed[1]]   # - 2\n",
    "    min_vals[:] = dist[seed[0]]              # - 3\n",
    "    upper_bound = min_vals.max()             # - 4\n",
    "    for n in seed:                           # - 5\n",
    "        np.min(np.array([min_vals, dist[n]]), axis=0, initial=upper_bound, where=np.logical_not(selected), out=min_vals)\n",
    "    # --- Loop argmax minimum ---\n",
    "    for n in range(n_seed, n_result):\n",
    "        sup_index = ma.array(min_vals, mask=selected).argmax()  # - 1\n",
    "        result[n] = sup_index                                   # - 2\n",
    "        v_dist[n - 1] = min_vals[sup_index]                     # - 3\n",
    "        selected[sup_index] = True                              # - 4     # | 5\n",
    "        np.min(np.array([min_vals, dist[sup_index]]), axis=0, initial=upper_bound, where=np.logical_not(selected), out=min_vals)\n",
    "    return result, v_dist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "执行上述代码后，就可以得到筛选出来的 4990 个样本 $\\boldsymbol{r}$ 序号，以及前 4889 个样本的采样距离 $\\boldsymbol{v}$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:06.217371Z",
     "iopub.status.busy": "2021-01-31T13:14:06.216605Z",
     "iopub.status.idle": "2021-01-31T13:14:06.833313Z",
     "shell.execute_reply": "2021-01-31T13:14:06.833716Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 628 ms, sys: 0 ns, total: 628 ms\n",
      "Wall time: 610 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(array([ 450, 4092, 3342, ..., 1696, 4495, 4400]),\n",
       " array([1963.36697 , 1819.286499, 1770.680538, ...,  959.175021,  956.828118,    0.      ]))"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "ks_sampling_core(dist, seed, n_result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## C 语言代码"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在 C 语言代码中只进行默认种子筛选与 Kennard-Stone 采样过程。我们暂且相信 numpy 可以很好地处理距离矩阵计算的过程。\n",
    "\n",
    "其编写是通过 ctype 的 Python/C binding 实现的。C 文件为 {download}`ks_cpp.c`，其 Python 接口文件为 {download}`KS_Sampling.py`。调用方式可以与上面相似："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:06.836698Z",
     "iopub.status.busy": "2021-01-31T13:14:06.836300Z",
     "iopub.status.idle": "2021-01-31T13:14:06.853589Z",
     "shell.execute_reply": "2021-01-31T13:14:06.853994Z"
    }
   },
   "outputs": [],
   "source": [
    "from KS_Sampling import ks_sampling_core_cpp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:06.857044Z",
     "iopub.status.busy": "2021-01-31T13:14:06.856599Z",
     "iopub.status.idle": "2021-01-31T13:14:06.916685Z",
     "shell.execute_reply": "2021-01-31T13:14:06.916321Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 158 ms, sys: 2.97 ms, total: 161 ms\n",
      "Wall time: 56.3 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(array([ 450, 4092, 3342, ..., 1696, 4495, 4400]),\n",
       " array([1963.366943, 1819.286499, 1770.680542, ...,  959.175049,  956.828125,    0.      ]))"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "ks_sampling_core_cpp(dist, seed, n_result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "如果不希望提供种子文件，且希望 C 程序自己找到两个距离最远的样本作为种子，那么通过下述方式调用："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:06.919543Z",
     "iopub.status.busy": "2021-01-31T13:14:06.919083Z",
     "iopub.status.idle": "2021-01-31T13:14:06.989368Z",
     "shell.execute_reply": "2021-01-31T13:14:06.988996Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 250 ms, sys: 11.4 ms, total: 261 ms\n",
      "Wall time: 65.5 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(array([4092,  450, 3342, ..., 1696, 4495, 4400]),\n",
       " array([1963.366943, 1819.286499, 1770.680542, ...,  959.175049,  956.828125,    0.      ]))"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "ks_sampling_core_cpp(dist, n_result=n_result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "如果想抽取所有样本，查看样本的抽样顺序，那么只需要提供距离矩阵即可："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:06.991453Z",
     "iopub.status.busy": "2021-01-31T13:14:06.991046Z",
     "iopub.status.idle": "2021-01-31T13:14:07.060941Z",
     "shell.execute_reply": "2021-01-31T13:14:07.061341Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 254 ms, sys: 8.13 ms, total: 262 ms\n",
      "Wall time: 65.6 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(array([4092,  450, 3342, ..., 1768, 1419, 3485]),\n",
       " array([1963.366943, 1819.286499, 1770.680542, ...,  921.688293,  910.627808,    0.      ]))"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "ks_sampling_core_cpp(dist)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在具体的 C 语言实现中，考虑了以下的优化代码因素：\n",
    "\n",
    "- 尽可能使用并行。在选取最大值、更新最小距离向量时，都使用 OpenMP 并行。并行的方式基本是默认。\n",
    "\n",
    "- 避免大量的指针地址计算。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 完整的 Kennard-Stone 采样函数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在引入上面定义过的函数 `get_dist`, `ks_sampling_core`, `ks_sampling_core_cpp` 的情况下，可以定义下述最终的采样函数；它的必选输入是样本数据，与其他类似程序的函数签名比较接近。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:07.063351Z",
     "iopub.status.busy": "2021-01-31T13:14:07.062925Z",
     "iopub.status.idle": "2021-01-31T13:14:07.067214Z",
     "shell.execute_reply": "2021-01-31T13:14:07.067610Z"
    }
   },
   "outputs": [],
   "source": [
    "def ks_sampling(X, seed=None, n_result=None, get_dist=get_dist, backend=\"Python\"):\n",
    "    \"\"\"\n",
    "    ks_sampling(X, seed=None, n_result=None, backend=\"Python\")\n",
    "    \n",
    "    KS Sampling Program\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    \n",
    "    X: np.ndarray, shape: (n_sample, n_feature)\n",
    "        Original data, need to be generated by user.\n",
    "        \n",
    "    seed: np.ndarray or list or None, shape: (n_seed, ), optional\n",
    "        Initial selected seed.\n",
    "        If set as `None`, the program will find the two samples\n",
    "        which have largest distance as seed.\n",
    "        \n",
    "    n_result: int or None, optional\n",
    "        Number of samples that should be selected.\n",
    "        If set as `None`, `n_sample` will be used instead, i.e.\n",
    "        selectet all data.\n",
    "    \n",
    "    get_dist: function\n",
    "        A function `get_dist(X)` that will read original data, and\n",
    "        return distance.\n",
    "        Default Implementation is Euclidean distance.\n",
    "    \n",
    "    backend: str, \"Python\" or \"C\"\n",
    "        Specify Kennard-Stone sampling function backend in Python\n",
    "        language or C language.\n",
    "    \"\"\"\n",
    "    X = np.asarray(X, dtype=np.float32)\n",
    "    if n_result is None:\n",
    "        n_result = X.shape[0]\n",
    "    dist = get_dist(X)\n",
    "    if backend == \"Python\":\n",
    "        if seed is None or len(seed) == 0:\n",
    "            seed = np.unravel_index(np.argmax(dist), dist.shape)\n",
    "        return ks_sampling_core(dist, seed, n_result)\n",
    "    elif backend == \"C\":\n",
    "        return ks_sampling_core_cpp(dist, seed, n_result)\n",
    "    else:\n",
    "        raise NotImplemented(\"Other backends are not implemented!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "如果希望在含有 5000 个数据的样本 `X` 中，依 Kennard-Stone 算法挑选 4990 个样本，可以用下述函数给出结果："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:07.069424Z",
     "iopub.status.busy": "2021-01-31T13:14:07.069017Z",
     "iopub.status.idle": "2021-01-31T13:14:07.812894Z",
     "shell.execute_reply": "2021-01-31T13:14:07.812542Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.32 s, sys: 47.6 ms, total: 1.37 s\n",
      "Wall time: 738 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(array([4092,  450, 3342, ..., 1696, 4495, 4400]),\n",
       " array([1963.367065, 1819.286499, 1770.680542, ...,  959.174927,  956.828064,    0.      ]))"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "ks_sampling(X, n_result=4990)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "更为快速的方法是使用 C 语言作为 Kennard-Stone 算法引擎："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:07.815914Z",
     "iopub.status.busy": "2021-01-31T13:14:07.815468Z",
     "iopub.status.idle": "2021-01-31T13:14:07.969356Z",
     "shell.execute_reply": "2021-01-31T13:14:07.969839Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 530 ms, sys: 44.5 ms, total: 574 ms\n",
      "Wall time: 151 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(array([4092,  450, 3342, ..., 1696, 4495, 4400]),\n",
       " array([1963.367065, 1819.286499, 1770.680542, ...,  959.174927,  956.828064,    0.      ]))"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "ks_sampling(X, n_result=4990, backend=\"C\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "如果希望前两个被采样样本序号是 (345, 456)，那么就需要需要指定种子大小 (采样顺序会严格地全部输出，但此时不一定会输出第一个样本距离)："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:07.971711Z",
     "iopub.status.busy": "2021-01-31T13:14:07.971309Z",
     "iopub.status.idle": "2021-01-31T13:14:08.113447Z",
     "shell.execute_reply": "2021-01-31T13:14:08.113843Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 514 ms, sys: 35.7 ms, total: 550 ms\n",
      "Wall time: 138 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(array([ 345,  456,  450, ..., 1696, 4495, 4400]),\n",
       " array([1388.464722, 1649.251587, 1633.396362, ...,  959.174927,  956.828064,    0.      ]))"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "ks_sampling(X, seed=[345, 456], n_result=4990, backend=\"C\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最后，如果希望计算的是 L1 范数距离而非欧氏距离 (即 L2 范数)，那么就更改 `get_dist` 函数即可："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:08.115700Z",
     "iopub.status.busy": "2021-01-31T13:14:08.115297Z",
     "iopub.status.idle": "2021-01-31T13:14:11.400109Z",
     "shell.execute_reply": "2021-01-31T13:14:11.398427Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.metrics import pairwise_distances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:11.404452Z",
     "iopub.status.busy": "2021-01-31T13:14:11.404050Z",
     "iopub.status.idle": "2021-01-31T13:14:12.611697Z",
     "shell.execute_reply": "2021-01-31T13:14:12.612030Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.28 s, sys: 36 ms, total: 1.31 s\n",
      "Wall time: 1.2 s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(array([ 345,  456,  450, ...,  999, 2095, 2046]),\n",
       " array([11158.767578, 13465.605469, 13155.375   , ...,  7544.62793 ,  7521.103516,     0.      ]))"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "ks_sampling(X, seed=[345, 456], n_result=4990, get_dist=lambda X: pairwise_distances(X, metric=\"l1\"), backend=\"C\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "之所以上述计算过程耗时长，是因为 L1 范数距离的耗时很大。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 代码效率比较"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下面是一些代码效率的比较，以示程序的高效程度与能力范围。\n",
    "\n",
    "这些代码是在无图形界面的 Linux 下运行，CPU 为 Intel-i5 7300HQ，内存为 16 GB，默认在可以 4 核并行时尽可能并行。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 距离矩阵生成 (5000×100 数据)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们这里比较三种做法。\n",
    "\n",
    "- `get_dist` 是上面的程序实现方法；\n",
    "\n",
    "- `get_dist2` 使用了更为直观的做法，即直接生成 $d_{ij} = \\sqrt{t_i - 2 y_{ij} + t_j}$；\n",
    "\n",
    "- `euclidean_distances` 是最为常用的生成欧氏距离的库函数之一。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:12.614714Z",
     "iopub.status.busy": "2021-01-31T13:14:12.614304Z",
     "iopub.status.idle": "2021-01-31T13:14:12.615750Z",
     "shell.execute_reply": "2021-01-31T13:14:12.616148Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.metrics import pairwise_distances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:12.618968Z",
     "iopub.status.busy": "2021-01-31T13:14:12.618568Z",
     "iopub.status.idle": "2021-01-31T13:14:12.619998Z",
     "shell.execute_reply": "2021-01-31T13:14:12.620399Z"
    }
   },
   "outputs": [],
   "source": [
    "def get_dist2(X):\n",
    "    Y = X @ X.T\n",
    "    t = Y.diagonal()\n",
    "    dist = np.sqrt(t[:, None] - 2 * Y + t[:, None])\n",
    "    return dist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`get_dist` 的耗时最少："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:12.623451Z",
     "iopub.status.busy": "2021-01-31T13:14:12.623049Z",
     "iopub.status.idle": "2021-01-31T13:14:16.006298Z",
     "shell.execute_reply": "2021-01-31T13:14:16.006718Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "161 ms ± 149 µs per loop (mean ± std. dev. of 7 runs, 3 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit -r 7 -n 3\n",
    "get_dist(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`get_dist2` 作为直观解法，耗时多一些，且内存消耗更大："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:16.008530Z",
     "iopub.status.busy": "2021-01-31T13:14:16.008131Z",
     "iopub.status.idle": "2021-01-31T13:14:20.719367Z",
     "shell.execute_reply": "2021-01-31T13:14:20.718938Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "224 ms ± 63.2 µs per loop (mean ± std. dev. of 7 runs, 3 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit -r 7 -n 3\n",
    "get_dist2(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`pairwise_distances` 库函数相对较慢，但仍然支持并行运算，效率相对还算可观。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:20.722757Z",
     "iopub.status.busy": "2021-01-31T13:14:20.722287Z",
     "iopub.status.idle": "2021-01-31T13:14:26.677504Z",
     "shell.execute_reply": "2021-01-31T13:14:26.677082Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "283 ms ± 5.63 ms per loop (mean ± std. dev. of 7 runs, 3 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit -r 7 -n 3\n",
    "pairwise_distances(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### KS 采样程序 (500×10 数据)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下面我们会考察一些具体的采样程序。我们会对网络上可以找到的开源代码做速度上的比较与测评，并简单验证我们 $O(n_\\mathrm{sample} n_\\mathrm{result})$ 程序的正确性。\n",
    "\n",
    "由于大多数开源的代码是 $O(N^3)$ 计算量级，且编程细节存在差异，因此只能在短时间内处理非常小的样本。我们先拿 500 样本的数据做验证。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:26.680460Z",
     "iopub.status.busy": "2021-01-31T13:14:26.680026Z",
     "iopub.status.idle": "2021-01-31T13:14:26.682707Z",
     "shell.execute_reply": "2021-01-31T13:14:26.682273Z"
    }
   },
   "outputs": [],
   "source": [
    "n_sample  = 500\n",
    "n_feature = 10\n",
    "X = np.random.randn(n_sample, n_feature)\n",
    "X *= 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下述代码是用 Python 程序，完整采样 500 个数据："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:26.685589Z",
     "iopub.status.busy": "2021-01-31T13:14:26.685171Z",
     "iopub.status.idle": "2021-01-31T13:14:26.713660Z",
     "shell.execute_reply": "2021-01-31T13:14:26.713058Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 74.8 ms, sys: 264 µs, total: 75.1 ms\n",
      "Wall time: 25 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "result_my_python = ks_sampling(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下述代码是用 C 程序，完整采样 500 个数据："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:26.716836Z",
     "iopub.status.busy": "2021-01-31T13:14:26.716276Z",
     "iopub.status.idle": "2021-01-31T13:14:26.749679Z",
     "shell.execute_reply": "2021-01-31T13:14:26.750174Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 93.1 ms, sys: 0 ns, total: 93.1 ms\n",
      "Wall time: 31 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "result_my_cpp = ks_sampling(X, backend=\"C\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "通过 `np.allclose` 可以两种代码的验证采样顺序是否一致。上面的 C 语言代码并没有更快，是因为小样本下 C 语言代码体现不出速度优势。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:26.753599Z",
     "iopub.status.busy": "2021-01-31T13:14:26.753107Z",
     "iopub.status.idle": "2021-01-31T13:14:26.756549Z",
     "shell.execute_reply": "2021-01-31T13:14:26.757032Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.allclose(result_my_python[0][2:], result_my_cpp[0][2:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下面找了比较常见的三个开源实现的纯 Python 代码并做比对；用于比对的程序放在 {download}`KS_Sampling_Others.py` 中：\n",
    "\n",
    "- https://hxhc.github.io/post/kennardstone-spxy/\n",
    "\n",
    "- https://github.com/karoka/Kennard-Stone-Algorithm/blob/master/kenStone.py\n",
    "\n",
    "- https://github.com/XiaqiongFan/PC-CCA/blob/master/PC-CCA.py\n",
    "\n",
    "第一个程序的实现效率较快，但尚没有达到我们上文定义的 `ks_sampling` 函数。在大数据量下，这种差异会更明显。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:26.759999Z",
     "iopub.status.busy": "2021-01-31T13:14:26.759505Z",
     "iopub.status.idle": "2021-01-31T13:14:26.787793Z",
     "shell.execute_reply": "2021-01-31T13:14:26.788366Z"
    }
   },
   "outputs": [],
   "source": [
    "from KS_Sampling_Others import ks_from_hxhc, ks_from_karoka, ks_from_XiaqiongFan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:26.791320Z",
     "iopub.status.busy": "2021-01-31T13:14:26.790493Z",
     "iopub.status.idle": "2021-01-31T13:14:27.070771Z",
     "shell.execute_reply": "2021-01-31T13:14:27.070324Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 651 ms, sys: 13.7 ms, total: 664 ms\n",
      "Wall time: 275 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "result_hxhc = ks_from_hxhc(X, test_size=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:27.074173Z",
     "iopub.status.busy": "2021-01-31T13:14:27.073693Z",
     "iopub.status.idle": "2021-01-31T13:14:27.076572Z",
     "shell.execute_reply": "2021-01-31T13:14:27.076153Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.allclose(result_hxhc[0][2:], result_my_cpp[0][2:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:27.079337Z",
     "iopub.status.busy": "2021-01-31T13:14:27.078932Z",
     "iopub.status.idle": "2021-01-31T13:14:33.457635Z",
     "shell.execute_reply": "2021-01-31T13:14:33.457196Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 6.91 s, sys: 71.9 ms, total: 6.98 s\n",
      "Wall time: 6.38 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "result_karoka = ks_from_karoka(X, n_sample)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:33.460811Z",
     "iopub.status.busy": "2021-01-31T13:14:33.460419Z",
     "iopub.status.idle": "2021-01-31T13:14:33.464141Z",
     "shell.execute_reply": "2021-01-31T13:14:33.463721Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.allclose(result_karoka[2:], result_my_cpp[0][2:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:33.678611Z",
     "iopub.status.busy": "2021-01-31T13:14:33.644735Z",
     "iopub.status.idle": "2021-01-31T13:14:42.668964Z",
     "shell.execute_reply": "2021-01-31T13:14:42.668541Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 9.26 s, sys: 114 ms, total: 9.38 s\n",
      "Wall time: 9.2 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "result_XiaqiongFan = ks_from_XiaqiongFan(X, n_sample)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:42.672376Z",
     "iopub.status.busy": "2021-01-31T13:14:42.671911Z",
     "iopub.status.idle": "2021-01-31T13:14:42.674842Z",
     "shell.execute_reply": "2021-01-31T13:14:42.674378Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.allclose(result_XiaqiongFan[0][2:], result_my_cpp[0][2:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### KS 采样程序 (2000×100 数据)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在选用一个适中的数据大小。在当前大小下，我们的程序仍然可以快速进行计算："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:42.677620Z",
     "iopub.status.busy": "2021-01-31T13:14:42.677230Z",
     "iopub.status.idle": "2021-01-31T13:14:42.683667Z",
     "shell.execute_reply": "2021-01-31T13:14:42.683994Z"
    }
   },
   "outputs": [],
   "source": [
    "n_sample  = 2000\n",
    "n_feature = 100\n",
    "X = np.random.randn(n_sample, n_feature)\n",
    "X *= 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:42.686812Z",
     "iopub.status.busy": "2021-01-31T13:14:42.686387Z",
     "iopub.status.idle": "2021-01-31T13:14:42.872496Z",
     "shell.execute_reply": "2021-01-31T13:14:42.872905Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 710 ms, sys: 15.8 ms, total: 726 ms\n",
      "Wall time: 184 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "result_my_python = ks_sampling(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:42.874673Z",
     "iopub.status.busy": "2021-01-31T13:14:42.874265Z",
     "iopub.status.idle": "2021-01-31T13:14:42.901850Z",
     "shell.execute_reply": "2021-01-31T13:14:42.902261Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 96.4 ms, sys: 0 ns, total: 96.4 ms\n",
      "Wall time: 24.2 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "result_my_cpp = ks_sampling(X, backend=\"C\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:14:42.904083Z",
     "iopub.status.busy": "2021-01-31T13:14:42.903681Z",
     "iopub.status.idle": "2021-01-31T13:15:02.824962Z",
     "shell.execute_reply": "2021-01-31T13:15:02.824523Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 19.4 s, sys: 1.1 s, total: 20.5 s\n",
      "Wall time: 19.9 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "result_hxhc = ks_from_hxhc(X, test_size=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:15:02.828125Z",
     "iopub.status.busy": "2021-01-31T13:15:02.827643Z",
     "iopub.status.idle": "2021-01-31T13:15:02.830879Z",
     "shell.execute_reply": "2021-01-31T13:15:02.830411Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.allclose(result_hxhc[0][2:], result_my_cpp[0][2:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### KS 采样程序 (30000×100 数据)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在相当大的数据集下，我们的程序可以确实地在较短时间内进行采样。对于个人电脑而言，由于 30000 个数据点所需要的距离矩阵已经相当庞大 (对于单浮点的数据，它占用大约 3.4 GB)，因此可以认为是个人计算机可以处理的数据量的极限。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:15:02.834479Z",
     "iopub.status.busy": "2021-01-31T13:15:02.833928Z",
     "iopub.status.idle": "2021-01-31T13:15:02.836846Z",
     "shell.execute_reply": "2021-01-31T13:15:02.836432Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.3527612686157227"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "30000**2 * 32/8 / 1024**3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在这个数据量下，我们的程序仍然能相当快地进行采样。对于高效且并行的 C 程序而言，其速度是非常快的；其最耗时的一步或许不是 Kennard-Stone 采样过程，而是距离矩阵的生成过程。而 Python 程序尽管较慢，但耗时仍然是可以接受的。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:15:02.839678Z",
     "iopub.status.busy": "2021-01-31T13:15:02.839292Z",
     "iopub.status.idle": "2021-01-31T13:15:02.914311Z",
     "shell.execute_reply": "2021-01-31T13:15:02.913876Z"
    }
   },
   "outputs": [],
   "source": [
    "n_sample  = 30000\n",
    "n_feature = 100\n",
    "X = np.random.randn(n_sample, n_feature)\n",
    "X = np.array(100 * X, dtype=np.float32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:15:02.917695Z",
     "iopub.status.busy": "2021-01-31T13:15:02.917298Z",
     "iopub.status.idle": "2021-01-31T13:17:21.231581Z",
     "shell.execute_reply": "2021-01-31T13:17:21.231155Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "19.8 s ± 10.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit -r 7 -n 1\n",
    "result_my_python = ks_sampling(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:17:21.234844Z",
     "iopub.status.busy": "2021-01-31T13:17:21.234403Z",
     "iopub.status.idle": "2021-01-31T13:17:57.690075Z",
     "shell.execute_reply": "2021-01-31T13:17:57.690480Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5.21 s ± 5.35 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit -r 7 -n 1\n",
    "result_my_cpp = ks_sampling(X, backend=\"C\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T13:17:57.692339Z",
     "iopub.status.busy": "2021-01-31T13:17:57.691930Z",
     "iopub.status.idle": "2021-01-31T13:18:18.359436Z",
     "shell.execute_reply": "2021-01-31T13:18:18.359011Z"
    },
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.95 s ± 2.24 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit -r 7 -n 1\n",
    "dist = get_dist(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[^Kennard-Stone.T.1969]: Kennard, R. W.; Stone, L. A. Computer Aided Design of Experiments. *Technometrics* **1969**, *11* (1), 137–148. doi: [10.1080/00401706.1969.10490666](https://doi.org/10.1080/00401706.1969.10490666)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
